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Abstract. 

In previous work with Scullard, we have defined a graph polynomial Pe(g, T) that 
gives access to the critical temperature T c of the (/-state Potts model defined on a 
general two-dimensional lattice C. It depends on a basis B , containing nxm unit cells 
of £, and the relevant root T c (n, m) of Ps(g,T) was observed to converge quickly to 
T c in the limit n, m oo. Moreover, in exactly solvable cases there is no finite-size 
dependence at all. 

In this paper we show how to reformulate this method as an eigenvalue problem 
within the periodic Temperley-Lieb algebra. This corresponds to taking m —> oo 
first, so that the bases B are semi-infinite cylinders of circumference n. The limit 
implies faster convergence in n, while maintaining the n-independence in exactly 
solvable cases. In this setup, T c (n) is determined by equating the largest eigenvalues 
of two topologically distinct sectors of the transfer matrix. Crucially, these two sectors 
determine the same critical exponent in the continuum limit, and the observed fast 
convergence is thus corroborated by results of conformal field theory. 

We obtain similar results for the dense and dilute phases of the 0(A) loop model, 
using now a transfer matrix within the dilute periodic Temperley-Lieb algebra. 

Compared with our previous study, the eigenvalue formulation allows us to double 
the size n for which T c (n ) can be obtained, using the same computational effort. We 
study in details three significant cases: (i) bond percolation on the kagorne lattice, up 
to n max = 14; (ii) site percolation on the square lattice, to ?i max = 21; and (iii) self¬ 
avoiding polygons on the square lattice, to ?i max = 19. Convergence properties of T c (n) 
and extrapolation schemes are studied in details for the first two cases. This leads to 
rather accurate values for the percolation thresholds: p c = 0.524404 999167439(4) 
for bond percolation on the kagome lattice, and p c = 0.592 746 050 79210(2) for site 
percolation on the square lattice. 


1. Introduction 

The question what makes a two-dimensional lattice model amenable to exact solution 
has attracted considerable attention within the field of statistical mechanics. Most, 
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but not all, solutions have been found by the technique of integrability, in which the 
commutativity of an infinite family of transfer matrices is ensured by requiring the 
Boltzmann weights to solve a set of cubic functional relations, known as the Yang- 
Baxter equations pL]. 

Recent work has focussed on a construction called discrete holomorphicity (DH), 
in which suitable correlation functions are required to satisfy a discrete version of the 
Cauchy-Riemann equations [2j. This leads to linear relations among the Boltzmann 
weights, that express the conservation of certain non-local currents in an associated 
quantised affine algebra M. The appropriate discretely holomorphic observables have 
been defined for several types of models, including Ising and Z N models [5], loop models 
of the Potts [6j and O (N) types [7], the chiral Potts model [8], and more exotic models 
involving multi-coloured loops [2]. 

In a series of papers with C.R. Scullard [TTH HT1 12] we have defined a topologically 
weighted graph polynomial for Potts and site percolation problems, having properties 
that are somewhat reminiscent of those found in the DH approach. This polynomial Pb 
depends on the Boltzmann weights of the degrees of freedom living within a “basis” B , 
by which we mean a small repeating part of the lattice. The main similarity of the graph 
polynomial approach with DH is, that when a set of Boltzmann weights corresponding 
to an exact solution is inserted, it produces a root of P B , independently of the size of 
B. In particular, when this size-independence is observed, it can be seen as heuristic 
evidence that we have found an exact solution. An important difference with DH is that 
the graph polynomial is defined as a partition function (albeit with some topological 
weighting of configurations), unlike the discretely holomorphic observables that take the 
form of correlation functions. 

The graph polynomial method is also practically useful when the model is not 
exactly solvable. Given a fixed set of physical coupling constants, let us denote by P B (T ) 
the evaluation of P B with Boltzmann weights corresponding to the temperature T. It 
is then observed maun ini m that the root T B —that is, a solution of P B (T B ) = 0— 
converges very quickly towards the critical temperature T c , upon increasing the size of B. 
This extends to models possessing several critical points [10], and even to inhomogeneous 
models with quenched bond disorder (spin glasses) [15]. 

The property of P B just mentioned can then be used as a numerical tool for 
determining T c very accurately. This was pursued extensively in [13] for the Potts 
model defined on all Archimedian lattices, their duals and their medials, as well as for 
site percolation on selected lattices (Archimedean and dual Archimedean lattices having 
only cubic and quartic vertices). 

The purpose of this article is to enhance the efficiency of this method, and to 
place it in a larger perspective by making contact with a number of existing theoretical 
constructions. To this end, we consider bases B consisting of n x rn unit cells of the lattice 
C on which the model is defined. According to mm , P B is defined by endowing B 
with doubly periodic boundary conditions and imposing a certain topological weighting 
of each configuration. The key idea in the present paper is then to take the m —» oo 
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limit first, so that the bases effectively become semi-infinite cylinders of circumference 
n. This transforms the criterion Pb(Tb) = 0 into an equality between eigenvalues of 
two topologically distinct sectors of the corresponding transfer matrix. These eigenvalue 
problems can then be solved—analytically for small bases, and numerically for larger 
ones—within the framework of the periodic Temperley-Lieb algebra [13] . 

This construction has several advantages. First, it makes the computation of 
Tb — T c (n ) numerically much more efficient—allowing basically for doubling the size n 
attainable, with respect to the previous approach ra — while maintaining the crucial 
feature that T c (ri) has no n-dependence at all when the model is exactly solvable. 
Second, it makes useful contact with both the transfer matrix formalism and with 
conformal field theory (CFT), allowing for a better understanding of the method. Third, 
it extends the applicability of the method beyond Potts and site percolation problems 
|13j . In particular, we shall show how to adopt it to O (N) loop models, in both the 
dense and dilute phases, in which case the underlying algebra is the dilute periodic 
Temperley-Lieb algebra. Fourth, since twice as many values T c (ri) are available for 
a given problem, the fini te-size scaling (FSS) behaviour can be studied much more 
carefully, and we can devise extrapolation schemes which are more efficient and reliable 
than the Bulirsch-Stoer acceleration of convergence employed in |13j . 

We illustrate all these aspects by applying the method to three significant unsolved 
problems, which in the past have each served as benchmarks within their respective 
category: 

(i) Bond percolation on the kagome lattice. This model has been the subject of a long¬ 
standing debate, since Wu’s ingenious 1979 conjecture for the percolation threshold 
p^ n = 0.524429 717 • • • [T5]. This conjecture was however proved incorrect, both by 
subsequent numerical work pa], among which p c = 0.52440499(2) [17] appears to 
be the most precise value to this date, and—maybe on a more fundamental level— 
by the observation m that Wu’s conjecture is exactly the outcome of the graph 
polynomial method with the smallest possible lxl basis (containing 6 edges). The 
previous graph polynomial method [l3j gave p c = 0.524404 999173(3), a result 
which we can now improve to 

p c = 0.524404 999167439(4). (1) 

This exemplifies the first and fourth points made above, revealing in particular that 
the error bar of [T3] was slightly underestimated. 

(ii) Site percolation on the square lattice. This renowned problem is unsolved essentially 
because the four-regular square-lattice hypergraph is not selfdual. The percolation 
threshold is in this case given by p c = 0.592 746 05(3) from numerical simulations 
S3. and by p c = 0.592 746 01(2) from the previous graph polynomial method [13]. 
We here improve this value to 

p c = 0.592 746 050 79210(2). (2) 

(iii) Self-avoiding polygons (SAP) on the square lattice. This is the polymer (N —* 0) 
limit of a dilute O (N) loop model on the square lattice, in which each vertex 
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can be visited at most once by the polygon. There is a fugacity 0 per monomer, 
and no bending rigidity. This model has been investigated extensively by exact 
enumeration techniques using CZQ. The best known critical monomer fugacity is 
z c = 0.379 052 277752(3) [21] . We have obtained for this problem values of z c (n) 
up to n max = 19. While the extrapolation of these data is compatible with—and 
more precise than—the result of m, we defer the discussion of extrapolations to a 
subsequent paper in which several numerical approaches to the SAP problem will 
be compared m 


The paper is organised as follows. In section [2] we review the definition of the graph 
polynomial for the g-state Potts model. The transfer matrix formalism, to be used 
extensively in this paper, is set up in section [3j Section [4] discusses the m —* oo limit 
that leads to the eigenvalue method for the Potts model. Using a few implementational 
tricks (see section [5]), we can then go on to determine the critical thresholds p c {n) for 
two selected percolation problems in section [6] The finite-size scaling behaviour of p c (n) 
is discussed in section [7] This leads to a powerful extrapolation scheme that provides 
the numerical values of p c = lim n ^. 0O p c (n) given in the abstract. Examining the relation 
of the eigenvalue method to conformal field theory, in section [8j enables us to generalise 
it to other representations (see section [9]) and other models. In particular, section 10 


sets up the method for the O (N) model, and discusses its relation with exactly solvable 


models. The application to the SAP problem is also provided there. Finally, section 11 
contains a few concluding remarks and perspectives for further investigations. 


2. Graph polynomial 

We first briefly review the definition of the graph polynomial for the case of the Potts 
model [[TO]. To set the scene for the eigenvalue method, we pay special attention to the 
ameliorations of computational complexity that were obtained in [12] 13] . 

Given a connected graph G = (V, E ) with vertex set V and edge set E, the partition 
function Z of the g-state Potts model [23] can be defined as m 

z = 27 «'V <A) , (3) 

AQE 

where |A| denotes the number of edges in the subset A, and k(A) is the number 
of connected components (including isolated vertices) in the subgraph Ga = (V,A). 
The temperature variable is denoted v = e K — 1, where K is the reduced interaction 
energy (including the inverse temperature) between adjacent g-component spins. In the 
representation (|3]) one can formally allow both g and v to take arbitrary real values. 
The special case of bond percolation is obtained by setting q = 1 and choosing the 
probability of an open bond as p — jf- } . 

The definition of the graph polynomial Pg(g, v) made initially in [10] was in terms of 
a deletion-contraction principle, whose validity it well-known for the partition function 
itself. The subsequent work [T2] provided an alternative definition that better reveals the 
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topological content of Ps(q,v). We henceforth suppose that G is an infinite, regular, 
two-dimensional lattice £. We define a basis B to be a finite subgraph of £ that 
produces all of £ upon application of an appropriate infinite set of translations that we 
call the embedding. The definition made in [42] is then in terms of conditioned partition 
functions, similar to (J3]) , that are defined on a graph G = (V, E) which is equal to the 
basis B: 

P B (q,v) = Z m - qZm . (4) 

Here Z 0D is the sum over edge subsets ACE such that all connected components 
(clusters) in the subgraph Ga = (V, A ) have trivial homotopy (i.e., are contractible to a 
point) upon endowing B with toroidal boundary conditions, whilst Z 2 d corresponds to 
clusters that wrap both periodic directions. Note that the terms in Z id, corresponding 
to clusters wrapping one but not the other periodic direction, do not appear in Q. 

From a computational point of view, the contraction-deletion algorithm described 
in [It)] has time and memory requirements that grow like 2^, where \E\ is the number of 
edges in B. In practice we shall be interested in bases that are iixm patterns of the least 
possible unit cell for B] a multitude of examples was given throughout [13]. For a regular 
lattice £ one has \E\ = kcnm, and for Archimedean lattices with the square embedding 
considered in [13], the proportionality constant kc ranges from 4 (square lattice) to 9 
(three-twelve and cross lattices); see Table 3 of [13] for details. In the transfer matrix 
algorithm of p2j the growth in time and memory is only like 4 2 ( n + m ) ; where 2 [n + m) 
corresponds to the perimeter (number of terminals) of B , that is, the number of vertices 
shared with translated copies of B within the embedding. Finally, Ref. [13] provided 
an improved transfer matrix within the periodic Temperley-Lieb algebra in which one 
of the periodic boundary conditions on B was imposed “on the fly”; this is shown in 
Figure 3 of m As a result, the exponential growth was reduced to 4 2mm ( ra > m ). 

In practical computations, one typically takes m = n, and so the computational 
effort is in nm. 256 n in [12], and 16 n in [13]. Accordingly, the maximum value of 
n that could be obtained for the Archimedean lattices was improved from n max = 2 in 
ma to n max = 4 in [13] . and further to n max = 7 in [13 ] . 

In this paper we show how to take the limit m —> oo, so that the basis B becomes 
a semi-infinite cylinder of circumference n unit cells of £. The roots of Pb(q, v) in that 
limit can then be computed by solving an eigenvalue problem in the periodic Temperley- 
Lieb algebra. The corresponding transfer matrix then acts on only n terminals, and time 
and memory requirements reduce to 4 n . Accordingly, the computations can now be 
taken to n max = 14 for the Potts model on the Archimedean lattices. We shall illustrate 
this below, by computing the bond percolation threshold on the kagorne lattice to high 
precision. Other values of q, and other lattices, are obviously also of interest, but a more 
systematic investigation will be reported separately [25] . 
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We shall refer to Q as the Fortuin-Kasteleyn (FK) representation of Pb(q, v). Each 
connected component in the subgraph Ga = (V, A) will be called an FK cluster. As in 
|13j . we shall need an equivalent formulation in terms of a loop model [26] defined on 
the medial lattice M(B). We now review the salient points leading to the definition of 
the transfer matrix in the loop representation, as well as the connectivity states that it 
acts on. 

The correspondence between FK clusters and loops can be depicted graphically as 
follows: 



(5) 


Here e = (ab) E E is an edge of G, and the left (resp. right) picture represents the 
situation where e E A (resp. e ^ A). The equivalent loops (shown in red colour) live on 
they bounce off the edge subset A and cut through its complement E\A. 

To turn this local equivalence into a global one, one uses the Euler relation for a 
planar graph to rewrite the partition function ^ in the loop representation as (26] 

z = d lVl/2 Y xlAl ‘ n iool , ( 6 ) 


ACE 


where x = v/y/q, and t(A) denotes the number of closed loops induced by the 
configuration A. The loop fugacity is ni oop = yfq- We shall use the parameters (q, v ) 
and (ni 00 p , x ) interchangingly. 

The loops on M(B) provide a representation of the Temperley-Lieb (TL) algebra. 
Supposing the direction of “time” propagation in ([5]) to be upwards, the left (resp. right) 
picture corresponds to the action of the identity operator I (resp. the TL generator E*) 
on the two adjacent strands, labelled i and i + 1 from left to right. We then have the 
relations m 


p2 _ p 

*“2 ^loopt-Z 5 

E, E ;; | E; E,. (7) 

E,E ; E,E ; for i - j > 1 . 


which may be proved graphically by gluing several diagrams in the form (J5]) on top of 
one another. 

The partition function ( 6b for a basis of size nxm can be computed within the TL 
algebra as shown in Figure [l The operator R* is here an element of the TL algebra built 
out of generators E ? with j E {i, i + 1, i + 2}, acting on connectivity states consisting of 
2 n strands, labelled 0,1,..., 2n — 1. The particular arrangement of Figure [l] is called a 
four-terminal representation of B |13j . Further details of this construction and explicit 
expressions for R,; for all Archimedean lattices can be found in [13]. 
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Figure 1 . Basis of size n x m in the loop representation. Terminals of the basis 
are shown as black circles, and periodic boundary conditions have been imposed 
horizontally. The auxiliary and quantum spaces, shown in red and blue colour 
respectively, sustain loops which are acted upon by an R^-matrix inside each grey 
square. 
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Figure 2. Three examples of connectivity states for n = 4. The numbers along the 
two time slices provide a canonical coding of the connectivity state. Loops are shown 
as red solid lines. The corresponding FK clusters live in the areas shaded in grey, 
whilst the dual FK clusters live in the white areas. 


A few examples of connectivity states are shown in Figure [2] The states are defined 
on two time slices (top and bottom), that describe the configuration of the system 
between vertical positions y = — 1/2 (bottom) and y = t — 1/2 (top), after any internal 
loop has been replaced by the factor ni oop , upon application of ([7]). The transfer matrix 
that propagates the system from “time” t to t + 1 is then given by the product of 
all R, operators within a row, where the periodic horizontal boundary conditions are 
implemented by tracing over the auxiliary spaces (shown in red colour in Figure [I]). 

Let us now be more specific about which variant of the TL algebra we actually 
need. The choice of periodic boundary conditions in the horizontal direction means 
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that the standard TL algebra [2?j (with generators E, for i = 0,1,..., 2 n — 2) must 
be made periodic. In the resulting periodic TL algebra there is an extra generator 
E 2n -i that acts across the periodic boundary condition (i.e., between strands 2 n — 1 
and 0). This algebra is however infinite-dimensional, and two modifications (algebra 
quotients) must be applied in order to make it finite-dimensional again. First, any loop 
of non-trivial homotopy (i.e., that winds around the periodic horizontal direction while 
being detached from the top and bottom time slices) must be replaced by a factor n win d. 
Second, loop segments connecting the two time slices are only considered according to 
which points they connect, and not how many turns they make around the periodic 
x-direction. The corresponding convention in Figure [2] is that among such connecting 
segments, the one that is leftmost in the top time slice (i.e., that carries the label ‘3’ in 
Figure [2]) is required not to cross the periodic direction. Note however that we still need 
to distinguish whether loop segments that connect a given time slice to itself crosses the 
periodic direction or not. With these modifications, the corresponding algebraic object 
is known as the augmented Jones- Temperley-Lieb algebra ; see section 3.1 of [28] for more 
details. 

In order to discuss further the states in Figure [2j we call arc a loop segment that 
connects two points within the same time slice, and string a loop segment that connects 
points on different time slices. The number of strings s = 2k is always even. When s > 0 
there are k FK clusters (and also k dual FK clusters) connecting the top and bottom 
time slices. A state can be turned into a pair of “reduced states” by cutting the s strings 
(if any) between the bottom and top time slices; each reduced state is then associated 
with only one time slice. Conversely, a pair of reduced states can be glued along the 
strings so as to form a “complete” state. This gluing can be done in k inequivalent 
ways, corresponding to cyclic rotates of the strings of one of the reduced states, in units 
of two (otherwise the distinction between FK clusters and dual FK clusters would fail 
to be respected). The s = 0 reduced states consist only of arcs and can be either closed 
(all FK clusters are bounded away from infinity) or open (at least one FK cluster is not 
bounded). For example, Figure [2 |d shows a state consisting of two open reduced states, 
whilst Figure^ depicts a pair of closed reduced states. Note that each reduced s = 0 
state can be conveniently coded in a binary convention where the code 1 (resp. 2) means 
a arc opening (resp. closing). 

A subtlety particular to the computation of Ps{q,v) arises because Z 1D does not 
appear in Q. We must therefore set n W i n d = 0. This implies that the s = 0 states can 
only be gluings of two open reduced states, or of two closed reduced states. In other 
words, mixed gluings are not allowed. We shall call such s = 0 states open or closed, 
respectively. 

Below, in section [4] we shall expose our main argument for transforming the 
computation of Prs(q, v) into an eigenvalue problem in the limit m —> 00 . On the 
level of the transfer matrix, this argument will imply an important simplification with 
respect to [18] . Namely, the eigenvalue problem can be solved by using only the top time 
slice, so that the transfer matrix acts only on reduced states. Moreover, only reduced 


Critical points of Potts and 0(N) models from eigenvalue identities 


9 


states without strings (s = 0) are needed—at least in the probabilistic regime v > 0 
considered here. There is an equal number of open and closed reduced states, namely 


1 

2 



~ 4 n 


( 8 ) 


of each. (In the intermediate stages of the computation, when n aux auxiliary spaces are 
open, simply replace n by n + n aux ). The gain of performance of the present method 
is largely due to the fact that ([8]) is much less than the number of “complete” states, 
which grows like ~ 16 n ; see Eq. (14) of [13]. 


4. Taking the m — > oo limit 

After these preliminaries, we now consider computing Pb(q, v) from Q for an n x m 
basis with finite n and m n; see Figure [T| In this limit, B has the geometry of a semi¬ 
infinite cylinder, which naturally suggests an interpretation in terms of the eigenvalues 
of a transfer matrix. 

When ordering the states according to a decreasing number of strings s, the transfer 
matrix T of section [3j with two time slices, has a lower block-triangular structure (the 
blocks being indexed by s), since under the time evolution the number of strings cannot 
increase. Its eigenvalues are therefore the union of the eigenvalues of each block on 
the diagonal, i.e., each eigenvalue can be characterised by the corresponding number of 
strings s. We denote these blocks by TW Moreover, the s = 0 block is a direct sum 
of two terms, T open and T c i oae< d, corresponding to a pair of open (resp. a pair of closed) 
reduced states. This is so precisely because contributions to Z id are excluded from Q, 
implying that loops winding around the cylinder carry the weight n wind = 0. So, as far 
as the eigenvalue problem is concerned, we can replace T by the direct sum 

n 

T = © T © T open © T closed . (9) 

k =1 

The modified transfer matrix T thus has the same spectrum as T, and moreover it 
conserves the quantum number s and, for s = 0, also the additional quantum number 
“open” and “closed”. Since T acts on the top time slice, and unlike T it cannot decrease 
s, it is unable to change the bottom reduced state. In other words, each of the direct 
summands T^ in ([ 9 J is in turn a direct sum of N s of identical blocks, with N s being the 
number of reduced states with s strings corresponding to the bottom time slice. Up to 
multiplicities of the eigenvalues, it therefore suffices to consider the action of (|9]) on the 
reduced states corresponding to the top time slice. In other words, the set of eigenvalues 
of T (that acts on “complete” states with two time slices) is the union of eigenvalues of 
T^ s \ T open and T closed , each restricted to act only on the set of reduced statesjf] 

We are interested in models with q > 0, and let us assume further that we are 
in the ferromagnetic regime, v > 0. All Boltzmann weights are therefore positive, and 


f A similar argument has been given in |29l 5JJ. 
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by the Perron-Frobenius theorem each summand in the decomposition ([9]) therefore 
has a unique, positive largest eigenvalue that we denote for s > 0, respectively 
A open and A c i ose d for s = 0. These eigenvalues obviously depend on the size n, and 
on the parameters (q,v). It follows from the cylinder geometry and the probabilistic 
assumption v > 0 that the eigenvalues are ordered 


-A-open? -^closed > A (2) > A |4) > • • • > A (2n) . 


( 10 ) 


Moreover, each of the terms in Q, Z 2 d and Z 0 d, must behave as ~ A m , where A is 
one of the above eigenvalues. By (10) the dominant contributions come from A open and 
Aciosed- Since Z 2 d (resp. Z 0 d) corresponds to an FK cluster (resp. a dual FK cluster) 
spanning the length of the cylinder, we must have 


^2D~(A open ) m , (11) 

(Aciosed)” 1 - (12) 


Note that these clusters will also wind around the circumference of the cylinder, with 
probability 1 in the limit m —>■ oo, since the number of strings is s = 0. 

It is obvious that A ope n and A c i OS ed are both greater than unity, so Z 2 d and Z 0 d grow 
exponentially with m. If (|4]) is to have a (positive, unique) zero as a function of v —as 
is indeed observed unmans] — there must exist some v c (n) > 0 so that A open = A c i ose d- 

We can prove this statement as follows. For r > 1 the dominant contribution to 
(J3|) will be A = E, and hence A open > A c i OS ed by direct computation. Conversely, for 
d< 1 the dominant contribution is A = 0, whence A operi < A c i OS ed. Since both terms in 
Q grow exponentially in m, the factor of q is unimportant, and the intermediate value 
theorem implies our main result 


Pr (diV) — 0 A open — Adosed i (13) 

valid for a basis B of size n x m, with n finite and m —> oo. 

We can therefore find the (unique) positive root, v > 0, of the graph polynomial for 
bases that are semi-infinite cylinders of circumference n by solving a simple eigenvalue 
problem over the s = 0 reduced states. Their number is given by ([8]) and grows as ~ 4 n , 
providing a substantial improvement over [13] in which Z 2 d and Z 0 d were computed by 
imposing complicated boundary conditions on a transfer matrix acting in the full set of 
states with dimension ~ 16 n . 


In the remainder of this section we first illustrate the main result (13) in a few 
simple cases. We end by discussing the special role of exactly solvable models. 


4-1. Square lattice with n = 1 

Consider the four-terminal representation of the square lattice with n — 1. There are 3 
reduced states which can be written ||, () and )(, where | denotes a string, ( is an arc 
opening, and ) is an arc closing. For the time being, we allow winding loops with weight 
riwi n d• It is easy to see that in the basis {||, (),)(} the transfer matrix reads T = T 2 T\ 
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T\ = 


T, = 


2x 

0 

0 



x 2 

^ wind 

2x + nioopX 2 

? 

(14) 

1 

2,X + Tl\ OQ1 p 

^wind 



2x 

0 

0 



1 

^wind 

2x + nioop 

1 

(15) 

x 2 

2x + nioop^" 

^wind*^ 



^-loop 

= yfq and x 

= vj-^fq. Setting n win d = 0 we obtain 



T = 


4x 2 0 0 

^loop T 4x (^ioop T 2x) 0 

x 4 (4 + ni 00 px) 0 x 2 (2 + nioopx) 2 


(16) 


which is indeed a lower block-triangular matrix, corresponding to the blocks s = 2, s = 0 
(closed), and s = 0 (open). Here, each of the blocks have dimension 1. The (dominant) 
eigenvalues read 


A (2) = 4 i 2 , Adced = ("loop + 2x) 2 , A open = x 2 (2+ n taa „x) 2 . 


(17) 


For n > 0 and x > 0 the ordering (10) is respected indeed 


To investigate the result (13) we remark that 

Aopen ^-closed ^■loop(^' l)(^hoop-£ T 4x T Ui 00 p) 


(18) 


is proportional to the graph polynomial Ps(q,v) = (v 2 — q)(q + 4v + v 2 ) for the 
n x m = 1 x 1 basis JTO] - In this case, where all blocks are one-dimensional, it is 
obvious—and can be verified by explicit computations—that the relevant (i.e., positive) 
roots of Ps(q, v) are independent of m. Note also that, despite of the assumption v > 0, 
this example actually correctly describes the phase diagram of the square-lattice Potts 
model both in the ferromagnetic (v > 0) and antiferromagnetic (v < 0) regimes, i.e., 
A< 2 > can simply be ignored. 


4-2. Kagome lattice with n — 1 

Consider next the kagome lattice with n — 1. Setting n W i nc j = 0 from the outset and 
considering only the s = 0 states {(),)(} we find that T is a diagonal 2x2 matrix with 
eigenvalues 

A dosed (("(loop T 2./() (xq CO p T 4/)] 00 p.x' T Gx T 2/),] oop .x' ), (19) 

Aopen = x 2 (2m oop + 12 a; + 13m oop x 2 + 6n 2 oop a; 3 + nf oop a; 4 ). (20) 

The difference A open — A closed is proportional to 

Ps(q, v) = v 6 + 6u 5 + 9n 4 — 2 qv 3 — 12 qv 2 — 6 q 2 v — q 3 , (21) 

which is the graph polynomial for the lxl basis. Back in 1979, Wu [15] conjectured 
this expression to be the exact critical manifold of the kagome-lattice Potts model, but 
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m 

T 

2 

4 

8 

16 

32 

oo 


Pc _ 

0.524429717521274793546879681534455071620567416578664793997510 

0.524406723188231819143234479992589885410333714096742273226669 

0.524406058417857416583229008103273638077164830301055000364284 

0.524406057896062955151905518860778390220248322088553687927465 

0.524406057896062634245378836787730760263849423348568828123454 

0.524406057896062634245378836666345666792028877197553627352494 

0.524406057896062634245378836666345666792028877197553609980248 


Table 1 . Bond percolation threshold p c on the kagome lattice using bases of size 2 x m 
for various m. 


our recent work PHJH21II3] definitively established that this is only an approximation 
corresponding to the smallest possible choice of the basis B. 

Once again the roots of Ps(q,v) for all the 1 x m bases are independent of m, 
because the blocks T open and T closed are one-dimensional. 


4-3. Kagome lattice with n — 2 


In Table [T| we show the estimates for the bond percolation threshold p c obtained as 
the relevant roots of Pg( l,v), with p = v/( 1 + v), for bases of size n x m with fixed 
n — 2 and varying m. The results for finite m were obtained from the algorithm of 


They are compared with the m — oo result obtained from (13) using the transfer matrix 
construction of the present paper. 

As expected, the results converge rapidly to the m = oo limit, the rate of 
convergence being exponential in m. We also note that the 2 x oo result 0.524406 057 • • • 
(i.e., the semi-infinite cylinder basis) is closer to the true percolation threshold p c = 
0.524404999 •• • |T3] than is the 2x2 result 0.524406 723 •• • (i.e., the square basis). 
Taking m —)■ oo in the 2 x m results is however far from “gaining one size”, since the 
3x3 result is 0.524405 172 ■ • These observations extend to arbitrary values of n. 


4-4- Exactly solvable cases 

It was shown in [13) that Ps{q,v) factorises over the integers for the three-terminal 
lattices (square, triangular and hexagonal), and that in the Ising case Pg(2, v) factorises 
for any lattice. In these factorisations, one or more “small” factors were observed to be 
independent of the size n x m of the basis, and their corresponding roots coincided 
with known exact solutions. In addition, a few sporadic cases were found in ra. 
mainly concerning Pg(0, v), where a size-independent factorisation occurred, and it was 
conjectured that these cases would be exactly solvable. 

Because of the m -independence of this factorisation result, it should still hold in the 
m —> oo limit. We therefore expect that A open — A closed will factorise in exactly solvable 
cases, spawning an n-independent factor whose roots provide the exact critical points. 
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We have already seen this happen in (18), and we have verified by explicit computations 
that this is indeed so also for higher n and for other exactly solvable models. 


5. Practical considerations 


To take the study to larger sizes n, we have implemented numerically the computation 
of the eigenvalues A open and A c i ose d that enter our main result (13). The transfer matrix 
T needs to be diagonalised in the sectors T open and T c i osed , and the first question to be 
settled is which is the most efficient technique for doing so. 

We have seen that T can be written as a product of R, operators, and these can in 
turn be written as product of the elementary operators 


H,: — I + xE i , Vj — x\ + Ej, (22) 

that add respectively a horizontal (or “space-like”) and a vertical (or “time-like”) edge 
to the lattice ra- The operators H, an V* are very sparse, with at most two non-zero 
entries per column, so the computation of W 2 — Tuj\ , where W\ is a vector of dimension 
Q, can be done with time and memory requirements which are proportional to that 
dimension. This calls immediately for iterative diagonalisation techniques |31j . 

The method of choice within this category is the Arnoldi method. However, we shall 
need to compute the eigenvalues to high numerical precision (cf. Table [Tj) , and we do not 
know of an implementation of the Arnoldi method which is compatible with arbitrary 
precision libraries. Moreover, we need only the largest eigenvalue in each of the sectors 
T open and T c i ose d, so a very simple method should be sufficient for our purposes. We 
have therefore used the most naive scheme, the so-called power method, in which the 
operations W 2 = Tw\ followed by w\ = w 2 /\\w 2 \| are iterated until | w 2 |j and w\ have 
converged to the largest eigenvalue of T and its corresponding eigenvector, respectively. 
To obtain convergence of the eigenvalue to 40-digit numerical precision, it turned out 
necessary to perform several hundreds of iterations, in particular for large sizes n. 

The action of Temperley-Lieb generators on the reduced states was described in [13] . 
Since we have s = 0 there are significant simplifications. The insertion and removal of 
the n aux = 2 auxiliary spaces (made necessary by the four-terminal representation of 
Figure [Tj) are handled exactly as in pT3] • In our implementation the states are stored in 
a hash table, since we want to take full advantage of the fact that for some problems 
(like, for instance, site percolation on the square lattice) the number of states needed 
can be even less than (J8[) . 

Another major simplification of the eigenvalue method is that no complicated 
topological considerations—such as those made in section 3.7 of [|13j— will be required 
in order to distinguish contributions to Z 2 d and Z 0D in Q. To choose the sector, it 
suffices to start the iterative scheme with an initial vector v equal to one of the reduced 
states in the ‘open’ or ‘closed’ sector, respectively. 

Considerations about efficiency are not limited to choosing the optimal numerical 
scheme for computing the eigenvalues. Once we can evaluate the function f(y ) = 
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A op en — A c i ose d for some value of the temperature variable v, we need also an efficient 
means of adjusting v to its critical value v 0 = v c (n ) satisfying f(v 0 ) = 0. Bracketing 
methods for finding zeros of a continuous function are numerically very stable, but rather 
slow. If we allow ourselves to compute derivatives, we can use instead the Newton- 
Raphson method, and more generally with fc’th order derivatives we can employ the 
fc’tli order Householder method. The higher-order methods will in principle converge 
faster, but in practice the computation of high-order derivatives is numerically unstable, 
so some compromise must be found. 

In practice we have found that the best result is provided by the second-order 
Householder method, with derivatives being computed by the symmetric difference 
method. Suppose we perform the computations in d-digit arithmetics (in practice we 
have taken d = 40). Let us set e = 10 -d/,2 _ Given some estimate v close to v 0 , one 
Householder iteration proceeds as follows. Make three evaluations of f(v), 

9o = f(v-e), 9i = f(v), 92 = f(v + e), (23) 


and apply the following formulae for the finite-difference derivatives: 

9-2 ~ 9o r (92 ~ 9i 


fo — 9u f i — 


2e 


/2 = 


9i ~ 9o \ -i 


(24) 


We stress here that to achieve numerical stability, the order of operations has to be 
carefully respected when computing f 2 . Finally, the next approximation to r; 0 is given 
by the second-order Householder formula 

fofi 


= v 


(/1) 2 - i/o/2 ' 

Extensive tests of this method shows that it has the following nice properties: 


(25) 


(i) After a few iterations, v converges to v 0 to full d-digit precision. 

(ii) If v is chosen sufficiently close to Vo, the number of correct digits will double in 
each iteration. 


The main computational effort obviously goes into the largest sizes n, and it is 
important in those cases to provide the best possible starting value 7j in i t for v. Thanks 
to the scaling theory developed in section [FJ we have been able to predict u init so that 
l^init — r>o| < HT 15 or better. This means that we can attain our d = 40 digit goal in 
just two Householder iterations. In a few cases we have contented ourselves with just a 
single Householder iteration and a final precision of at least 30 digits. 


6. Percolation thresholds on selected lattices 

In this section we apply the eigenvalue method to two prominent sample problems which 
have been extensively studied in the past: site percolation on the square lattice, and 
bond percolation on the kagome lattice. 

We stress that it would be straightforward to study also the Potts model for other 
values of q ^ 1, or to switch to any of the other lattices treated in j!3j . To this end, it 
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suffices to change the R; matrix to any of the explicit expressions given in ra. which is 
a matter of changing just a few lines of code. We shall however defer these extensions to 
a future study [25], and for the time being take the computations for the two problems 
mentioned as far as possible. 

6.1. Site percolation on the square lattice 

The R-matrix for site percolation on the square lattice is given by Eq. (76) in [13]: 

R = Ej + 2 E ; + uE l+ i, (26) 

where we recall that q — 1 and the probability of an open bond is p = This 

R-matrix contains only two out of fourteen possible terms, so thanks to the use of 
hashing techniques only a subset of the reduced states will be used in the diagonalisation 
procedure (see sections 3.6 and 7.0 of H3) for more details). 

The site thresholds were found on nxn square bases up to n max = 11 in [13j. Using 
the eigenvalue method we have obtained the thresholds on n x oo bases up to n max = 21. 
These results are shown in Table [2j A comparison with Table 52 of [13] shows that the 
10 x oo result is very close to the old 11 x 11 results, so having taken the m —> oo limit 
can be said, roughly speaking, to have “gained us one size” in this case. 

It is obvious from Table [2] that there is agreement with the existing results for p c 
on at least the first five digits. More digits can however be obtained by extrapolating 
the data, and this will be discussed in section [7] 

6.2. Kagome lattice (3, 6, 3, 6) 

The R-matrix for the Potts model on the kagome lattice is given by Eq. (29) of [13]: 

R Hi + iVj-|_ 2 VjE i+ iV i+ 2VjHj + i, (27) 

where the elementary operators H, and V, were defined in ( ]22| . We set q = 1 to obtain 
the corresponding bond percolation problem. 

The bond thresholds were found on nxn square bases up to n max = 7 in [13]. 
From the eigenvalue method we have computed the thresholds on n x oo bases up to 
u-max = 14. Those results are shown in Table [3j The value p c (2) was already presented 
in Table [H 

It is immediately visible that the data in Table [3] converge faster than those in 
Table [2j The value of p c (n ) with n = 14 agrees with p c to eleven digits. Looking down 
the table we can also see that the central value for p c given in Ref. ra is a bit too high, 
and that the final extrapolated result is likely to be slightly below its lower error bound. 
We shall discuss this extrapolation in section [7] 

7. Extrapolations 

We shall now discuss the extrapolation of the data in Tables [2]-[3] in view of obtaining 
final values of the thresholds p c which are as precise as possible. 
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n 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 

Ref. pi] 
Ref. 1 13] 


Pc(n) _ 

0.5000000000000000000000000000000000000000 

0.5651977173836393964375280132470308160984 

0.5888806999178529980514426957517049337221 

0.5914171708531384817988341017359231779642 

0.5922358232050266776468513240523872931777 

0.5925073562056416791039647019136652231541 

0.5926196333998949001725078647635478154618 

0.5926727605746273159803396143033787878155 

0.5927006240698093405431688044620515383831 

0.5927163956307984449936472582379354676976 

0.5927258706594202658220083530779420439346 

0.5927318424282431711880689641407018936896 

0.5927357579756109329667635225705210166147 

0.5927384119896431219655807545828528768826 

0.5927402625774169696977289423157479414628 

0.5927415848750128489249007208681102997016 

0.5927425500481430481174634605214209833281 

0.5927432678876343617903095425293033143864 

0.5927438107312915517933469085441350226515 

0.5927442273849199618209333613184304919328 

0.592744551481371482002735520463 

0.59274605 (3) 

0.59274601 (2) 


Table 2. Site percolation thresholds p c (n ) on the square lattice, as computed from 
n x oo bases, and two previous results for p c . 


7.1. Site percolation on the square lattice 


A reasonable Ansatz for the asymptotic behaviour of p c (n), motivated by general 
principles of fini te-size scaling (FSS), is that of a series of power-law corrections, 


Pc{n) = p c + 



k=1 


(28) 


with 0 < Ai < A 2 < • • •. To determine Ai we first form 8p c (n ) = p c (n) — p ci where p c 
is taken either as the existing best value pn Eg, or from a preliminary fit to the data 
of Table [2} We then consider the sequence 


A / x = log (8pc(n)) - log (5p c (n - 1)) 
1 U log(n) — log(n — 1) 


( 29 ) 


which by (28) should converge to the first FSS exponent Ai. 
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n 

Pc(n) 

1 

0.5244297175212747935468796815344550716205 

2 

0.5244060578960626342453788366663456667920 

3 

0.5244050922187183914064917102789956045159 

4 

0.5244050138823434506779249332748912013263 

5 

0.5244050026660985339974686380437997371046 

6 

0.5244050002521386411660652383853120094089 

7 

0.5244049995708026048576486896416033403853 

8 

0.5244049993387487061840419066777093506317 

9 

0.5244049992479802098067018389586796534330 

10 

0.5244049992084754512628552771194320896813 

11 

0.5244049991897559735118013090108129282307 

12 

0.5244049991802484437799696383468246717858 

13 

0.5244049991751328450538203030184876090453 

14 

0.524404999172242908087780703763 

Ref. p3] 

0.52440499 (2) 

Ref. [TBj 

0.524404999173 (3) 


Table 3. Bond percolation threshold p c on the kagome lattice, as computed from 
nxoo bases, and two previous results for p c . 



1/n 

Figure 3. Determination of the first FSS exponent Ai for site percolation on the 
square lattice. 
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Figure 4. Determination of the second FSS exponent A 2 for site percolation on the 
square lattice. 


In Figure |3] we show Ai(n) as a function of l/n. The data are very well fitted by 
a polynomial in 1 /n 2 , and allowing for some freedom on the degree of the polynomial 
and the number of small-n points to be excluded from the fit, we arrive at the result 

A 1 = 4.000 1(2). (30) 

This agrees well with the value w = 4.03±0.01 reported in section 7.2 of ra. which was 
found as the optimal choice for the FSS exponent entering the Bulirsch-Stoer algorithm. 
It appears inevitable to admit that Ai = 4 exactly. 

We next form the sequence 

n 4 p c (n) — (n — 1 ) A p c {n — 1) 


p A (n) = 


n 4 — (n — l) 4 


(31) 


in which the leading FSS term A\/n 4 has been subtracted off (28). From this we form 
a sequence 

log (5p A (n)) - log (Sp 4 (n - 1)) 


A 2 (n) = 


(32) 


log(n) — log(n — 1 ) 

which should now converge to the second FSS exponent A 2 . We plot A 2 (n) against l/n 
in Figure [I], along with a polynomial fit in l/n 2 . This yields 

A 2 = 6.00(1). (33) 

and we henceforth admit that A 2 = 6 exactly. 

It is now obvious how to continue. In the next round we subtract Ai/n 4 ± A 2 /n 6 
from (28) and seek to determine A 3 from the residue. Going through the same steps as 
above we find A 3 = 8.0(5), and we conjecture that A 3 = 8 . 

The determinations of A 1; A 2 and A 3 provide compelling evidence that A k = 
2 (k ± 1) for any k. The precise FSS form then reads 

OO 4 

A k 


Pc{n ) = Pc + 


k =1 


77,2(fc+l) 


(34) 
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1.5 

1 

D- U 

5° 

3 ? 

0.5 

0 

Figure 5. The estimators Pg"° (in arbitrary units) plotted against no- 



We now show how to use this form to obtain a very precise extrapolation for the 
percolation threshold p c . 

Let n max denote the largest size for which we have been able to compute p c (n). 
In the present case we have n max = 21, as seen in Table [2} We first form a series of 
estimators Pm,l hi which the scaling form (34) is truncated at the l/n M term, and in 


which the data p c (n) is used up to a maximum size of n = L. In other words, we find 
the unique solution of the linear system 

Am/ 2-1 

“ 1-JT T * * * - 

with n = L + 1 — M/2,..., L — 1, L. Second, for a fixed M, we form another series 


'A 1 A 2 

PM,L + [ -4 + “6 + 


= Pc(n), 


(35) 


of estimators p'/ff } by fitting pm,l to the residual dependence predicted by (34), but 
eliminating from the fit the first n$ possible values of L. That is, we find the unique 
solution of the linear system 


rf" o) + 




+ 


Bo 


+ ••• + 


By, 


a—no—1—M/2 


Pm,l ■ 


(36) 


This is a fit on n r 


^M+4 ^2(?T m ax n 0 1) 

no — M/2 different values of L ranging from 1 + M/2 + no up to 


n„ 


If we eliminate too few data points (i.e., take no too small) when forming the 
estimators p^ , the result will be mediocre because it depends too much on the smallest 


sizes for which the FSS form (34) is dubious. On the other hand, if we eliminate too 


many data points (i.e., take no too large) the result will again deteriorate because the 
fit has too few terms. We would expect an optimum in between these extremes. 

In Figure [ 5 ] we show the variation of p^ with no in the case M = 6. For reasons 
of clarity, we actually plot the quantity 10 11 {p^ ~ Pc), where p c is our final value 
for the percolation threshold, but the units of the ordinate in the plot should really 
been thought of as arbitrary, since we have not determined p c yet. We see that there 
is an extremum (minimum) at some intermediate value of n 0l in agreement with the 
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M 

Estimate 

4 

0.592746050791752 

6 

0.592746050792111 

8 

0.592746050792085 

10 

0.592746050792096 

12 

0.592746050792125 

14 

0.592746050792165 

16 

0.592746050792226 


Table 4. Estimates \ {'p'fff + p'^f ) for the site percolation threshold on the square 
lattice p c . 


above qualitative argument. Repeating the plot for other values of M (not shown), it is 
observed that the minimum becomes more shallow upon increasing M, at least up to a 
certain point beyond which the quality of the plot deteriorates due to a lack of points. 

For any value of M in a reasonable range (namely M = 4, 6,..., 16) the minimum 
is found to occur for no — 8 or uq — 9. The arithmetic mean of those two values 
of pffl can thus be taken as a precise estimate for p c . We show these mean values in 
Table [4] to 15 significant digits. They are seen to depend only very weakly on M in some 
intermediate range (say M — 6,8,10,12) from which we can extract our final value and 
error bar for the percolation threshold: 

p c = 0.592 746 050 79210(2). (37) 


7.2. Bond percolation on the kagome lattice 


For bond percolation on the kagome lattice we again start by considering an FSS Ansatz 


of the form (28). It is immediately clear from the data that the leading FSS term does 


not correspond to the exponent Ai = 4, as was found for site percolation on the square 
lattice. To obtain a unified notation we therefore set A x = 0 in the kagome case, and 
call the leading FSS correction A 2 /n A2 . 

Figure [6] shows estimates A 2 (n) for this leading FSS exponent—extracted from the 
data p c (n) just as in ( [29] )— , plotted against 1 jn. The accompanying fit is a second-order 
polynomial in 1/n, and this and similar fits lead to the value 

A 2 = 6.00(5). (38) 

With data up to only n max = 7, it was concluded in section 4.3 of ra that w ~ 6.36 
was a suitable exponent for the Bulirsch-Stoer extrapolation. It is clear from Figure [6] 
how this conclusion could be reached, since effectively only the the last six data points 
or so (we have now n max = 14) start bending upwards. Although this determination 


of Ao is somewhat less accurate than (33), we are confident in concluding that A 2 = 6 
exactly. 
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Figure 6. Determination of the leading FSS exponent A 2 for bond percolation on the 
kagome lattice. 


Comparing the precisions of A 2 and A 3 obtained in section 7.1 it is clear that there 


is little hope of obtaining a convincing determination of A 3 in the present case. However, 
the FSS exponents A fc are not only a property of these lattice models of percolation, but 
are also expected to characterise the field theory describing their continuum limit. There 
is overwhelming evidence throughout the literature that both models are described, 
in the continuum limit, by the same conformal field theory, which can be derived by 
standard Coulomb gas arguments [32]. So there are good reasons to believe that we 
should have A^ = 2 (k + 1) also in the present case. Different lattice realisations will 


however give different values of the non-universal amplitudes in (28). It is quite 


possible that the three-fold rotational symmetry of the kagome lattice (which replaces 
the four-fold symmetry of the square lattice) has the effect of setting A\ = 0. This also 
matches observations made in 


We therefore proceed with the analysis as in section 7.1, using in particular the 


scaling form (34) with A\ = 0. Going through the same steps as before we arrive at the 


final value for the percolation threshold 

p c = 0.524 404 999 167439(4). 


(39) 


8. Connection to conformal field theory 


The fact that the leading exponent in the scaling form (34) takes a high value, Ai = 4, 


is responsible for the fast convergence of p c (n) towards p c and the ensuing precise 


determinations (37) and (39). We now examine how this fast convergence can be linked 


to considerations about the continuum limit. 

The free energy per unit area fo{n) of a conformally invariant system defined on a 
semi-infinite cylinder of circumference n scales like 

/o(n) = /o(°°)-j^ + o(i 


C 2 ) • 


(40) 











Critical points of Potts and 0(N) models from eigenvalue identities 


22 


where c is the central charge of the corresponding conformal held theory (CFT) and 
/o(oo) is the bulk free energy. This can be related to the largest eigenvalue A 0 of the 
transfer matrix for a corresponding lattice model as fo(n) = — ^logAo, where £ is a 
geometrical factor that depends on the lattice (£ = 1 for the square lattice) and ensures 
the correct normalisation per unit area in the lattice model. 

Similarly, the free energy fi(n) of excited states (i — 1, 2,...) has the scaling [85] 


27ray 


fi(n) ~ fo(n) = + o (n 2 ) , 


(41) 


where Xi is the corresponding scaling dimension (critical exponent). The smallest 
excitation fi{n) of the CFT corresponding to percolation is related to the magnetic 
exponent x m , so we have X\ = x m . The values c = 0 and x m = A are of course known 
ra. but they are not important for the following argument. 

The important point is that the two transfer matrix sectors, T open resp. T c i ose d, 
considered in section [4] correspond to excitations in which an FK cluster (resp. a dual 
FK cluster) is required to propagate along the semi-infinite cylinder. Define now the 
corresponding free energies per unit area 


fopenip) log A open , 


/closed (jl) — lo§ A c ] ose d 


(42) 


In the continuum limit there is no difference between whether the propagating cluster 
is an FK cluster or a dual FK cluster. Therefore / ope n(^) and ./dosed ( n ) both determine 
the same critical exponent, namely x mi and they both scale like fi(n) in (41). It follows 
that the difference 


/open (p) /closed(^) O (fl ) (4*8) 

vanishes fast as n —> oo, right at the critical point p = p c . 

This is a suggestive argument, but it does not quite explain the convergence 
properties of the eigenvalue method. What we have observed in sections[4]and[7]is, that if 
we define a pseudo-critical point p c (n) as the value of p for which / op en(^) — /closed (n) = 0, 
then 


Pc(n) - p c = O (n 4 ) , (44) 

and moreover the corrections appear to be 0 (n -6 ), 0(n~ 8 ), and so on. 

It is clear that more work would be required to establish whether (43) can be 
shown—obviously using more ingredients—to actually imply (44). But one thing that 
has become clear is, that the eigenvalue method owes its success to the fact that / open (n) 
and /dosed (n) are two different ways of determining the same critical exponent. This 


will be exploited further in section 10 


9. Spin representation 

It is of interest to review the definition (j4j) of the graph polynomial Pb(q, v) when 
q G N. In that case the Potts model can be defined directly in terms of g-component 
spins, instead of the FK clusters that we have considered this far. 
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Let again the basis B consist of n x m unit cells of the lattice C. Define Z fi v as the 
partition function on B with doubly periodic boundary conditions that are twisted by 
/i = 0,1,..., q — 1 (resp. v — 0,1,..., q — 1) in the horizontal (resp. vertical) direction. 
By this we mean that the values of a pair of nearest-neighbour spins, cr* and a 3 , that 
are on opposite sides of the horizontal (resp. vertical) periodic boundary condition are 
considered identical if eq — aj = p mod q (resp. cq — crj = v mod q), and different 
otherwise. 

To relate the partition functions Z jiv in the spin representation to those in the FK- 
representation (Z 0D , Z 1D and Z 2 d) we first notice that untwisted boundary conditions 
are simply doubly periodic, whence 


Z oo — Z 0 d + Zid + Z 2 d • (45) 

Consider next the quantity 'Y^ I1V Z^ V . Configurations in Z 0D contribute to all q 2 
terms in this sum, whereas those in Z 2 r> can only contribute to one term, namely Zoo- 
Finally, configurations in Z id are such that all clusters that are non-homotopic to a 
point have the same topology, i.e., they have the same winding numbers ( n x ,n y ) with 
respect to the horizontal and vertical periodic boundary conditions [36J. These winding 
numbers are defined up to a global sign change, (— n xi —n y ) = ( n x ,n y ), and if they are 
both non-zero they must satisfy 


gcd (n x ,n y ) = 1 . 

Now, for a configuration in Z 1D to contribute to Z yu we should have 
n x p + n y v = 0 mod q . 


(46) 


(47) 


Thanks to the constraint (46), this equation has precisely q solutions for the labels 
(/i, v). Summarising, we have proved that 
q -1 9-1 

EE Z l:l , — q 2 Zot) + qZiv + Z 2 D ■ 


(48) 


/i=0 u=0 


Combining (45) and (48) we obtain 

q-l q -1 , , 

2»-iEE Z thV — f 1-J [Z 2 d — qZon) , 


q 


(49) 

fi=0 v=0 

so the quantity on the left-hand side is proportional to the graph polynomial Pe(g, v ) by 
Q. This was already shown in the appendix of P3] , by using a more involved argument 
of duality transformations. 

We now consider the m — > oo limit of (49) in order to obtain an eigenvalue criterion 
in the spin representation which is equivalent to Q. The asymptotic behaviour, as 
m —>■ oo, of the partition functions reads 

z„,„ =c„„(A„r+ dll (Apr+■ ■ ■, (bo) 

where the eigenvalues depend only on /i, but the coefficients can depend on both 
twist labels. We have ordered the eigenvalues in decreasing order: > A? > 
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Moreover, the dominant eigenvalues in each sector decrease when the twist increases, 
Aq > Ai > ..., where obviously A p = A g _^. Moreover, some of the inequalities might 
not be sharp when q takes particular values. 

Going back to the left-hand side of (49), we see that the dominant contribution 
(A 0 ) m cancels out between the two terms. The next-leading contributions come from 
A« and Ai. These must cancel out in order for Pniq, v) to vanish: 

P B (q,v) = 0 Ag 1} = Ai, (51) 

valid for finite n, in the limit m — > oo. This is the spin-representation version of our 
main result (13). 

The two eigenvalues involved have a very precise meaning. The leading eigenvalue 
Ao in the untwisted (periodic) sector corresponds to an eigenvector which is invariant 
under a global permutation of the spin, Oi — » per* with p G S g . The next-leading 
eigenvalue A^ transforms non-trivially under such a transformation: it picks up a non¬ 
trivial r/tli root of unity. For instance, when q — 2, Aq 1 " 1 is the largest eigenvalue 
corresponding to an eigenvector which is odd under spin reversal. It is well-known that 
the free-energy gap between Aq 1 ^ and Aq determines the magnetic exponent x m by (41). 


On the other hand, Ai is the largest eigenvalue in the twisted sector p — 1, in which 
spin labels are shifted cyclically by one unit when one crosses the periodic boundary 
condition. It is equally well-known that its free-energy gap with respect to the ground 
state A 0 determines the same exponent x m . It follows that the criterion (51) can be 
discussed in exactly the same terms as in section [8j 


10. O(N) loop model 

The graph polynomial method, and its eigenvalue version pursued in the present paper, 
can be seen as a tentative to generalise the notion of self-duality to situations, where 
duality is not an exact symmetry. In the Potts model, a duality transformation exists 
directly on the lattice, and interchanging the lattice and the dual lattice amounts 
to shifting cyclically the sites in the loop representation by one lattice unit. This 
transformation provides a bijection between the states in the open and closed sectors, 
as discussed in section |3l 

It is clearly of interest to formulate the graph polynomial method also for other 
models, and in particular for the O (N) loop model [57115%] . Although the two models are 
in the same universality class—to be more precise, the dense phase of the O (N) model 
has the same central charge and a closely related, albeit not identical, operator content 
as the critical g-state Potts model with N = ^Jq = m oop [32]—their lattice definitions 
present subtle differences, which were remarked early on |36j . The O(N) model does not 
possess a duality transformation on the lattice, and it treats parity issues in a different 
way than the loop formulation of the Potts model. In particular, in the periodic transfer 
matrix formalism, the 0(A r ) model can be defined on any number of strands, whereas 
in the Potts model the number of strands needs to be even (e.g., there were 2 n strands 
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in Figure [Tj) . Also, on a semi-infinite cylinder with free boundary conditions at the 
infinities, the number of non-contractiblc (winding) loops in the Potts model must be 
even, but in the O(N) model this number can have any parity. 

These subtleties have some important consequences in the continuum limit, as can 
be seen by examining the operator content of the models in detail. Arguably the most 
important difference is that the energy operators of the two models do not coincide, as 
can be seen from a detailed Coulomb gas (CG) analysis [SUES, 32J- The same analysis 
also reveals that the O (N) model possesses two kinds of involutions that could be viewed 
as duality symmetries. The first involution exchanges the dense and dilute theories 
corresponding to the same central charge, by replacing the CG coupling constant g by 
1/g. The second involution originates from the study of modular invariant partition 
functions [HJj and amount to exchanging the the role of electric and magnetic charges 
in the CG. 


10.1. Eigenvalue method 


The observations made in section [8] give us an important hint about how to obtain 
an eigenvalue criterion for the O(N) loop model, similar to (13) and (51) for the 
Potts model. It seems that we should try to identify one same critical exponent that 
arises in the continuum limit of two topologically distinct sectors of the transfer matrix. 
Moreover, the introduction to section 10 provides the clue that the two sectors should 
differ by their charge content (i.e., electric versus magnetic) in the Coulomb gas analysis. 

The algebraic framework of the 0(N) model is that of the dilute TL algebra. The 
precise setup that we shall need is that of the dilute augmented Jones-Temperley-Lieb 
algebra, whose description is as in section [3| except that we should now allow for dilution, 
in the sense that some vertices and edges are not covered by loops. To make this 
statement precise, we first recall that the Potts model defined on a planar graph G is 
equivalent [26] to a completely packed loop model defined on the 4-regular medial graph 
A4(G), with the loops around a vertex of M.(G) being in any of the two states ([5]). In 
the case of the 0(N) model, the loops are defined directly on the chosen graph G, which 
hence needs not be 4-regular [37], but to keep things simple we shall concentrate on the 
O(N) model defined on the square lattice [38]. The two states ([5]) are then replaced by 
the following nine states of the loops around a vertex 



We define R, as the sum over those nine diagrams, each one being weighed by the 
corresponding Boltzmann weight pi as shown. Integrable choices of R; exist [38] and 
will be discussed in section |10.2[ but for the moment we are interested in the general— 
and not necessarily critical—case where pi take arbitrary values. 

In this dilute model the TL generators E; are defined as before [see (J5])] , except 
that Ej can only be applied when the strands i and i + 1 are occupied by loop segments, 
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m 



H-1-1-1- *■ X 

0 1 ••• n — 1 


Figure 7. Transfer matrix construction in the loop representation for the O(N) 
model. Periodic boundary conditions have been imposed horizontally. The auxiliary 
and quantum spaces, shown in red and blue colour respectively, sustain loops which 
are acted upon by an R;-matrix inside each grey square. 


corresponding to the ninth diagram in (52). The algebraic relations (J7]) are unchanged, 
but must be supplemented by the requirement that each of the nine generators (52) can 
act only when the loop strands have the correct occupancy. This is graphically clear, 
but slightly cumbersome to write down in algebraic terms 05 - Note that the weight of 
a closed loop is now denoted N = n\ oop . 

The transfer matrix T is again defined as the product over R,;, followed by a trace 
over the auxiliary space, as shown in Figure [TJ As usual in the the algebraic approach 
to integrable systems, the direction of time propagation is upwards in diagrams, such 
as (52), defining the action of algebra generators, whereas time flows to the North-East 
when auxiliary spaces are present, ft follows that the diagrams (52) must be rotated 
45° in the clockwise direction before being placed inside the gray squares in Figure [7j 

The states on which T acts can still be drawn as in Figure [2j except that empty 
vertices are now possible (and can be represented by the special code ‘0’). 

There are two types of operators in the CG. A magnetic operator m s (r) inserts a 
topological defect, so that s oriented loop strands are created in a small neighbourhood 
around the point r. The two-point function (m s (ri)m_ s (r 2 )) corresponds, in the 
cylinder geometry where ?r and r 2 reside at either extremity, to imposing the 
propagation of s strings in the transfer matrix setup of section [3j In particular, the 
largest eigenvalue A( s ) of [cf. (19j)] determines the critical exponent x s of the operator 


m s , via (41). Note that s can have any parity in the O(N) model, unlike the Potts case 
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where s = 2k must be even. If the loop weight is parameterised as 
N = — 2 cos(7rg), 

where g is the CG coupling constant, then 

1 , il-9? 


x s = -as 
8 y 


2 g 


(53) 


(54) 


The other type of CG operator is the electric operator e e (r), also known as a vertex 
operator. Its key property is that the two-point function (e e (r 1 )e_ e (r 2 )) amounts, in the 
cylinder geometry and in the s = 0 sector, to setting the weight of each non-contractible 
loop to 

fVwind = 2cos(7re). (55) 

The particular choice e = e 0 = 1 — g corresponds to the usual situation iV wind = N, 
which provides the ground state of the model. The corresponding charge eo is called 
the background electric charge. For general e, the critical exponent with respect to the 
ground state reads [35] 

-(1-S) 2 


X P = 


2 g 


(56) 


The dense (resp. dilute) phase of the O (N) model corresponds to the regime 
0 < g < 1 (resp. 1 < g < 2). The results for the dense phase also apply to the 
critical Potts model, by setting N = ^fq. In particular, setting jV win d = 0 amounts 
to forbidding winding loops, so—by the reasoning of sections [3] and [8]— we obtain the 
magnetic exponent as x m = X\/ 2 - 

Motivated by the introductory remarks in this subsection, we now consider the 
lowest magnetic excitation mi, corresponding to having one string propagate along the 


cylinder, with exponent x\ given by (54). The electric exponent x e in (56) can be made 


to take the same value upon making a particular choice of the charge e: 

y n 

x P — x\ = - -— = 0 <=> 


e g n 

x\ — -= 0 

2 g 8 


,9 

6 — zlz — 


(57) 


This is equivalent to choosing 

N wind = ±y/2^N. (58) 

More generally, we would get x e = x s for e = ±sg/2, but taking the clue from the Potts 
result, we should focus on the closest equivalent of the magnetic (order parameter) 
operator in the Potts model, which is indeed mi in the O (IV) case [ 39 ]. 

Based on this argument, we define the eigenvalue method for the O (N) model as 
follows. For finite size n, find the value of the parameters p t (n) so that the largest 
eigenvalue in the s = 1 sector, coincides with the largest eigenvalue A in the s = 0 


sector with the particular choice (58): 
A (1) = A 


(59) 


with JVwind = ±a/2 — N. 

The sign ambiguity on the right-hand side will be resolved later. 

This proposed method succeeds or fails depending on whether it can deliver both 
features that distinguished the eigenvalue method for the Potts model: 
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(i) The values Pi(n) should be independent of n in exactly solvable cases. 

(ii) For noii-solvable cases, Pi(n) should converge “very fast” in n. 

This success criterion will be examined in details in the remainder of this section. 
But let us note for now one encouraging observation. The dense 0(1) model with 


vertices (52) is equivalent to a site percolation problem on the square lattice, with 


certain local interactions depending on p l . It is known that for a particular choice 
of pi , that corresponds to the integrable model |38j with an arbitrary inhomogeneous 
choice of spectral parameters, the ground state has a combinatorial nature that can be 
investigated 02 via the quantum Knizhnik-Zamolodchikov approach. The ground state 


in the s = 1 sector is also combinatorial, and in particular (59) holds true for any finite 
n, with the choice N wind — N — 1. 

10.2. Exactly solvable cases 

We now consider the integrable case [3B 1 ) 13j of the model (52) with weights 

sin(w) sin(3A — u ) 


pi(u) = 1 + 


P 2 (u) = p 3 (u) = 
Pa(u) = p 5 (u ) = 
Pe{u) = Pr{u) = 
Ps(u) = 


sin(2A) sin(3A) 
sin(3A — u) 


sin(3A) 
sin(w) 
sin(3A) 

sin(w) sin(3A — u) 
sin(2A) sin(3A) 
sin(2A — u ) sin(3A — u ) 


Ps(u) = - 


sin(2A) sin(3A) 
sin(u) sin(A — u) 


(60) 


sin(2A) sin(3A) 

The spectral parameter u governs the anisotropy of the interactions, and we have here 
written pi = Pi(u) for later convenience. The crossing parameter A is related to the loop 
weight via 

N = — 2 cos(4A). (61) 

We take arbitrary inhomogeneous spectral parameters, meaning that u = Uk for any 
vertex in the fc’th column of the lattice. 

For size n = 1 there is just one reduced state in either of the sectors s = 0 and 
s — 1. The one-dimensional transfer matrices T^ read 

T (0) = pi(wi) + N wind p 6 (u i), 

T (1) = p 7 (ui) +pa(ui) +p 9 (tti) • (62) 

Using trigonometic identities, the difference T (0) —T (1) is proportional to iV W md + 2 cos 2A, 


so we conclude that (59) is satisfied with 


= — 2 cos(2A). 


(63) 
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As A goes from 0 to tt, the loop weight N runs through the range [—2, 2] four 
times, corresponding to the four branches of critical behaviour discussed in [38] . The 
first two branches, A G [0, |] and A G [f , ^], correspond to the dilute and dense phase, 


respectively. The relation between (61) and (63) is such that the plus (resp. minus) sign 


in front of the square root in (59) should be taken for A G [|, (resp. A G [0, w] U [^, 7r]). 


In particular, the plus (resp. minus) sign should be taken for the dense (resp. dilute) 
phase of the O(N) model. 

We now seek confirmation of these results for size n = 2. In the sector T there 
are 3 reduced states which can be written oo, () and )(, where o denotes an empty site, 
and the remainder of the notation is as in section Id.ll In the sector T d) the 2 reduced 
states are |o and o|. With this ordering of the bases, the transfer matrices read 

pip'i + PqP&N p 3 p' 4 N + p4p' 3 N p 3 p' 4 N + p 4 p' 3 N 
P5p2 P7P7 + P9p' 8 N PsPs + P9pg + P9P 8 N 

P 2 P 5 Psp's + p9pg + PsPqN p 7 p' 7 + pSpgN 


7"(°) = 


(64) 


and 


j'W — 


P 7 P 1 + PsPq + P 9 P 6 
P'iPg + P4p' 5 


p2p3' + p b p' A 
P 1 P 7 + PePs + pep'g 


(65) 


where we have abbreviated piiu-f) = Pi, Pi(u- 2 ) = p\ and -/V wind = N. Inserting now 


(60), (61) and (63) we find that the two eigenvalues of T ^ coincide with two of the 


eigenvalues of T [) > , for arbitrary values of the parameters A, u\ and u 2 - 

We have similarly studied this model at size n = 3, in which case dim(T'h) = 7 
and dim(T^ 1 - > ) = 6. Remarkably, we found that with arbitrary inhomogeneous spectral 
parameters, u 2 and u 3 , and for arbitrary values of A, all 6 eigenvalues of T ( h were 
also eigenvalues of 

The dimensions of these transfer matrices are related to the Motzkin numbers. 
Define M(x) = (1 + x)(l — 3x) and consider the generating functions 

1 OO 

fo(x) = 7 == = 


y/Mfr) 


h(x) = 


= 2_^ a n x 

n =0 

2x 


Y. b > 


X 


( 66 ) 


71=1 


M (x) + (1 — x) a jM(x) 

Then a n = dim(T^) and b n = dim(T (1 ^) for a system of size n loop strands. We have 
a n > b n for n > 1. However, both numbers exhibit the same asymptotic behaviour for 
n » 1: 

- f—V 3 n . (67) 


“'/t ''ii 1 1 ~ 

2 \nnJ 

We conjecture that with inhomogeneous spectral parameters, all eigenvalues of T ^ are 
also eigenvalues of provided the weight of non-contractible loops is taken as in 


(63). If true, this would be very promising for finding a genuine graph polynomial for 


the O (N) model, i.e., one having properties similar to those of Pb(q, v) in the Potts case 
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[HH H2 US] for finite n x m bases, and not just in the m —> oo limit. We hope to report 
more on this soon. 


10.3. Approximation method 

We now investigate the second aspect of the eigenvalue method for the O (N) model, 
namely its usefulness as an approximation method for the critical points of non-solvable 
models. To this end, we apply it to the problem of self-avoiding polygons (SAP) on the 
square lattice, which is the N —y 0 limit of a loop model in which each occupied edge 
has the weight z. There is no bending rigidity, and the osculating vertices p 8 and pg are 


disallowed. The Boltzmann weights (52) are thus 


Pi = 1 , 

Pi = P3 = Pi = Pb — Pe — Pi — z , 

ps = P9 = 0 . (68) 

Note that we have suppressed the spectral parameter u, since this model is not 
integrable. 

This SAP model has been extensively studied by exact enumeration techniques 
[muni ed], and the critical monomer fugacity is known to very high precision iza 

z c = 0.379052 277752(3). (69) 

This value corresponds to the smallest z > 0 for which the generating function, as 
obtained by exact enumeration, exhibits a singularity. In the formulation of the problem 


in terms of a partition function, with the weights (68), this corresponds to selecting the 
dilute branch of the 0(N) model, and implies taking the minus sign in ( |59| ). We therefore 
set N = 0 and N wm( \ = — \/2. 

By diagonalising the transfer matrices T ^ and T W and proceeding as in section [HJ 
we have obtained values of z c (n) up to n max = 19 for the SAP problem. They are 
shown in Table [5] It is clear that these data exhibit the required fast convergence, and 
a detailed analysis—to be presented elsewhere m —reveals that the scaling is in fact 


compatible with (34). 


We should emphasise that the formulation in terms of a partition function allows us 


to study the problem (68) also for other values of N, and in the dense phase. The dilute 


N — 1 case corresponds to an Ising model on the square lattice in which configurations 

of alternating spins (H-1— or —I-h) around a lattice face have been disallowed, 

since we have set p 8 = p 9 = 0. This produces a value of the critical domain wall fugacity, 
z c = 0.421 326 ■ ■ ■, which is slightly lower than that of the standard Ising model, which 
reads ^ sing = (1 + v 7 ^)" 1 = 0.414214 ■ ■ ■. 

The dense phase of SAP (with IV = 0) does not appear to have been studied 
explicitly within the exact enumeration framework, but the corresponding z c is likely to 
manifest itself as a subdominant singularity of standard, dilute SAP. 
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n 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 
19 

Ref. gH 


z c (n) _ 

0.3832870437289217825415444959209990643484 

0.3800152822923947541103727449094743052839 

0.3793419092420152604076859124268482909456 

0.3791615386298805591124869699564102732536 

0.3791017465104568577033096312174651793134 

0.3790779263723816763857349117326710080035 

0.3790669419366251682820022783255996752011 

0.3790612863965732376129739341339159714858 

0.3790581237478262657302859193323348704028 

0.3790562392439348634338963536547147709970 

0.3790550583590770828697993099179842253186 

0.3790542873705249946097446478792002255473 

0.3790537664746062070854620937409756594548 

0.3790534041836437725305420784870138649786 

0.3790531458388626510867578645132848654379 

0.3790529575762840825464391257666224019613 

0.3790528177462476184521578790271596607432 

0.3790527121228867470 

0.379052277752 (3) 


Table 5. Critical fugacities z c (n ) of the SAP model on the square lattice, as computed 
from n x oo bases, and a previous result for z c . 


We finally note that when applying (59) to a non-solvablc model, the parameter 
z c (n) can indeed be tuned so that the leading eigenvalues of T® and T^ 1 ) coincide, but 
the remaining eigenvalues of T-- 1 will in general not be equal to eigenvalues of 7^°k 


11. Discussion 

In this paper we have transformed the graph polynomial method of uni hot into an 
eigenvalue method. This corresponds formally to taking the m —> oo limit of the nx m 
bases B that enter the definition of the graph polynomial Pg. The advantages of this 
reformulation are numerous and have been discussed in the Introduction. We can add 
to this list that, on a technical level, the eigenvalue method requires only the reduced 
states in the transfer matrix setup (see section [3]) and avoids the rather complicated 
topological considerations of [13j (see section 3.7 of that reference in particular), two 
facts that make the practical implementation of the method considerably easier. 

On a more fundamental level, the eigenvalue formulation has revealed that the 
method hinges on identifying two distinct topological sectors of the transfer matrix 
that lead to the determination of one same critical exponent. This has enabled us to 
extend the applicability of the method from the g-state Potts model—including bond 
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and site percolation problems—to encompass O (N) models in various phases, even in 
the presence of inhomogeneities. Also in the O (N) case we have demonstrated that the 
method is both capable of 

(i) detecting exact solvability (in the sense that integrable models lead to results 
independent of the size n), and of 

(ii) generating approximations to the critical parameters that converge rapidly in n. 


This first aspect poses a set of fundamental questions that should motivate future 
research. In particular, the possible link between exact factorisation in the graph 
polynomial method and discrete holomorphicity (alias conservation of non-local currents 
in quantised affine algebras 001), or related manifestations of exact solvability, remains 
to be elucidated. With the present extension from Potts to O (N) models, we have 
demonstrated that the graph polynomial method is likely to be as ubiquitous and 
versatile as discrete holomorphicity. The fact that the O (N) version of the method 
required us to impose a very particular value of A r wind in (59), and the ensuing massive 
eigenvalue coincidences between T ^ and T^°\ are strongly reminiscent of phenomena 
encountered in representation theory |4Tj, and this possible link should be examined as 
well. 


The second aspect has enabled us to study the finite-size scaling properties of the 
method in much more detail than m- In particular, we have developed powerful 
extrapolation schemes capable of determining the critical point to within 15-digit 
precision. In future work, we plan to extend these determinations to other lattices 
(following [13] in the Potts case), and to other values of the parameters q and N. 
From a more practical perspective, we are working on a parallel implementation of the 
algorithm which should make accessible larger n and lead to even higher precision ra. 

Finally, the extension to yet other models, such as Zn models and multi-coloured 
loops, previously considered from the discrete holomorphicity perspective 13 El, should 
also be investigated. The question of whether the present method applies to loop 
models with non-trivial boundary interactions jiol m U] provides another appealing 
perspective. 
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